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Abstract. 

We analyze three simple genetic circuits which involve transcriptional regulation and feed- 
back: the autorepressor, the switch and the repressilator, that consist of one, two and three 
genes, respectively. Such systems are commonly simulated using rate equations, that account 
for the concentrations of the mRNAs and proteins produced by these genes. Rate equations are 
suitable when the concentrations of the relevant molecules in a cell are large and fluctuations 
are negligible. However, when some of the proteins in the circuit appear in low copy numbers, 
fluctuations become important and the rate equations fail. In this case stochastic methods, such 
as direct numerical integration of the master equation or Monte Carlo simulations are required. 
Here we present deterministic and stochastic simulations of the autorepressor, the switch and 
the repressilator. We show that fluctuations give rise to quantitative and qualitative changes in 
the dynamics of these systems. In particular, we demonstrate a fluctuations-induced bistability 
in a variant of the genetic switch and and noisy oscillations obtained in the repressilator circuit. 

1. Introduction. The production of proteins in cells is regulated by networks of in- 
teracting genes. Some of these genes code for transcription factors, which are proteins 
that regulate the expression of other genes by binding to specific promoter sites on the 
DNA. Some of these proteins, called repressors, perform negative regulation, while others, 
called activators, perform positive regulation. Post-trancriptional regulation can occur by 
translational regulation, or by post-translational regulation such as protein-protein inter- 
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action, which may modify the function of these proteins. These networks of interacting 
genes are at the heart of all processes in cells: from the regulation of the cell cycle to the 
various stress responses. 

It turns out that genetic networks are typically sparse networks, namely most genes 
interact only with a small number of other genes. Also, the networks exhibit some degree 
of modularity, namely one can identify recurring modules, which are known as network 
motifs [18] , 

Modules found in nature are often hard to study and fully control. To overcome these 
limitations, synthetic networks can be constructed from well-studied components. Ex- 
amples of such components are the lactose, the lambda and the tetracycline systems in 
bacteria. The synthetic networks are designed to perform desired functions, determined 
by their architecture. They do not require the manipulation of the structure of proteins 
and other regulatory elements at the molecular level. The genes and promoters are of- 
ten inserted into plasmids rather than on the chromosome. Two important examples of 
synthetic circuits are the genetic toggle switch and the repressilator. 

The toggle switch consists of two genes, which negatively regulate each other expres- 
sion. The regulation is performed at the transcriptional level, namely each gene codes 
for a repressor protein that binds to the promoter of the other gene. Such system may 
exhibit bistability, namely two stable states, where in each state one of the proteins is 
dominant and the other is suppressed. A synthetic toggle switch was constructed in E. 
coli and the conditions for bistability were examined [B]. The switching between its two 
states was demonstrated using chemical and thermal induction. More recently, such cir- 
cuit was found to exist in a natural system in which two mutual repressors regulate the 
differentiation of myeloid progenitors into either macrophages or neutrophils 

The synthetic repressilator circuit, constructed in E. coli, was designed to exhibit 
oscillations, reminiscent of natural genetic oscillators such as the circadian rhythms |4J. 
The repressilator circuit was encoded on a plasmid that appears in a low copy number. 
The protein concentrations were measured vs. time in single cells. Oscillations with a 
period of about 150 minutes were found, namely extending over several division cycles. 
The oscillations were found to be noisy, typically maintaining phase coherence for times 
of the order of a single oscillation period. 

The dynamics of genetic networks is often simulated using rate equation models. These 
are sets of coupled ordinary differential equations, which account for the concentrations 
of the mRNAs and proteins in the network. In general, rate equations account for aver- 
age concentrations and ignore fluctuations. They are suitable for systems in which the 
concentration of interacting molecules are large and fluctuations are negligible. However, 
proteins in cells often appear in low copy numbers and may exhibit large fluctuations. 
Moreover, in case of transcriptional regulation, the expression of the regulated gene may 
be controlled by a single protein that binds to its promoter. This extends the notion of 
low copy numbers and the resulting fluctuations even to the case when there is a large 
number of free repressors of a certain type, since only one of them may bind to the pro- 
moter at any given time. Recent advances measurement techniques made it possible to 
measure the fluctuations in copy number of proteins in single cells [5J [331 US] ■ Measure- 
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ments of protein levels in single cells revealed distributions that depend on the topology 
of the regulatory network controlling the particular protein |19j . 

To account for the fluctuations, simulations of genetic networks should be done using 
stochastic methods such as the master equation [101 HH [TT1 [20l EU [22J, [23] or Monte 
Carlo simulations [9l I16j . The master equation provides the probability distribution 
of the concentrations of proteins in a population of cells. From this distribution one 
can calculate the average concentrations as well as the correlations and the rates of all 
processes. Monte Carlo simulations enable to follow the temporal variations in the level 
of gene expression and in the concentrations of proteins in a single cell. These results 
enable to extract the noise level and temporal correlation functions. They can also be 
used to characterize the dynamics of the system and determine whether it exhibits a 
single steady state, multi-stability or oscillatory behavior. 

In this paper we present deterministic analysis (using rate equations) and stochastic 
analysis (using direct numerical integration of the master equation and Monte Carlo 
simulations) of three simple genetic circuits: the autorepressor, the toggle switch and 
the repressilator (Fig. []}. We show that fluctuations give rise to both quantitative and 
qualitative effect in the dynamics of thse circuits. 

The paper is organized as follows. In Sec. 2 we consider the autorepressor and present 
the deterministic and stochastic methods used in the paper. In Sec. 3 we study the genetic 
switch and in Sec. 4 we analyze the repressillator circuit. The results are discussed and 
summarized in Sec. 5. 

2. The Autorepressor The autorepressor is the simplest genetic circuit that involves 
feedback. It consist of a single gene, denoted by a, that negatively regulates its own ex- 
pression. The gene transcribes into mRNAs that translate into A proteins. These proteins 
function as repressors. When an A protein binds to the a promoter site, it prevents the 
RNA polymerase from binding to the promoter and thus inhibits the transcription pro- 
cess. It turns out that genetic networks include a large number of autorepressor modules. 
This circuit is thus considered as a network motif [25) . One may speculate that it per- 
forms some crucial function in the cell. It was proposed that the role of the autorepressor 
is to speed up response times |25) or to reduce fluctuations 2 J . Below we analyze the au- 
torepressor using deterministic and stochastic methods. We utilize its simplicity in order 
to present the methodologies in details. 

2.1. Michaelis-Menten equations The dynamics of genetic networks are commonly de- 
scribed using the Michaelis-Menten equations. These equations describe the temporal 
variations in the concentrations of the relevant molecules in the cell. Here, we denote the 
concentration of protein A in a cell by [A] (by concentration we refer to the average copy 
number of A proteins per cell). The concentration of the corresponding mRNA is denoted 
by [m]. The Michaelis-Menten equations for the concentrations take the form 



M 



1 + k[A] n 
g p [m] - d p [A}. 




(1) 



[A] 
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Figure 1: The three genetic circuits considered in this paper: the autorepressor (a); the 
genetic switch (b); and the repressillator (c). The flat-headed arrows denote negative 
transcriptional regulation. 

The parameters d m and d p (sec -1 ) are the degradation rates of the mRNAs and proteins, 
respectively. The trascription rate is given by g m (sec -1 ). The rate in which each mRNA 
is translated into proteins is given by g p (sec -1 ). Since A proteins negatively regulate their 
own synthesis, the transcription rate is reduced by a factor of 1/(1 + fcL4] n ). This factor 
is called the Hill-function. In this expression, the parameter k quantifies the regulation 
strength (determined by the affinity between the repressor and the promoter site). The 
exponent n is called the Hill-coefficient. In general, Hill-function models can be derived 
from more complete rate equation models. In this case, n is expected to take only integer 
values. In fact, n represents the number of copies of the transcription factor, that are 
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required to be bound simultaneously to the promoter in order to perform the regulation. 
The case of n > 1 is often referred to as cooperative binding. In this paper we consider 
only integer values of n. However, similar models may also be used to fit empirical data. 
In this case, n is simply a fitting parameter which may take non-integer values. 

To simplify the analysis of genetic circuits, the mRNA level is often ignored and the 
transcription and translation processes are regarded as a single step of protein synthesis 
|25j . In this case, the effective production rate of proteins is given by g = g p g m /dm- The 
Michaclis-Menten equations are reduced to 

where d — d p . Ignoring the mRNA level is typically justified under steady state conditions. 
However, in the analysis of systems that are away from steady state due to external 
signals, or those that exhibit oscillations or large fluctuations the mRNA level should be 
included. The Michaelis-Menten equations for the autorepressor exhibit a single steady 
state solution for the concentration of A proteins. In case that n = 1, it is given by 

(3) [A\ - Zl±vgj§g . 

This solution is found to be stable for any choice of the parameters. 

2.2. Extended set of rate equations In the Michaelis-Menten equations the negative reg- 
ulation process is described by the Hill-function. This description is rather crude and 
incomplete. In order to model the regulation process in greater detail we present below 
a more complete rate equation model [12] • In this model, we account separably for the 
populations of free and bound proteins. The bound A proteins are denoted by r and their 
concentration is given by [r]. 

Consider an autorepressor gene a, encoded on the chromosome, which exhibits no 
cooperative binding. In this case the number of bound repressors is in the range < 
[r] < 1. The concentration, [r], can also be considered as the fraction of time in which the 
promoter is occupied by a bound repressor and the transcription process is suppressed. 
Therefore the transcription rate is reduced by a factor of (1 — [r]). Ignoring the mRNA 
level, the extended set of rate equations takes the form [12] 

[A] = g (i-[ r ])-d[A]-ao[A](l-[r]) + ai[r] 

(4) [r] = ao[A](l-H)- ai [r], 

where the parameter ao (sec -1 ) is the binding rate of the repressors to the promoter site 
and at (sec -1 ) is the their unbinding rate. In the limit in which the binding and unbinding 
processes are much faster than other processes in the system (namely ao, a% 3> d, g), these 
equations can be reduced to the Michaelis-Menten form. In this limit, the relaxation 
time of the concentration [r] is much shorter than other relaxation times in the system. 
Therefore, one can take the time derivative of [r] to zero, even if the system is away 
from steady state. This brings the rate equations to the Michaelis-Menten form [Eq. ^] 
with n = 1 and k = ao/a\. Therefore, Eqs. ^ have the same steady state solution 
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for the protein A as the Michaelis-Menten equation ([2|). However, the dynamics leading 
to steady state may differ between the two equations. Furthermore, the extended rate 
equation model exhibits more flexibility in the sense that it is much easier to insert 
additional features into Eqs. ((4|) than into Eq. ([2]). For example, it is possible to introduce 
degradation of bound repressors by adding the term — d r [r] to the equation for [r] in Eqs. 

©■ 

2.3. Stochastic Analysis Transcription factors and other proteins, as well as their mRNAs 
in a cell, often appear in low concentrations [9l [T71 [8] . In this case, fluctuations in the 
copy numbers of these molecules may play an important role in the dynamics of genetic 
networks. To obtain a better description of thse systems, one should take into account the 
discrete nature of the molecules rather than using continuous concentrations. Moreover, 
even in case that some transcription factor appears in a high concentration, the regulation 
is performed by a small number of copies that are bound to the promoter site. The 
fluctuations in the number of bound transcription factors give rise to large temporal 
variations in the transcription rate of the regulated gene. 

In order to account for fluctuations in genetic networks, stochastic methods are re- 
quired, such as the master equation or Monte Carlo simulations [U [10l [16j HQl H3] . The 
master equation is expressed in terms of the probability distribution P(Na, N r ). This is 
the probability for a cell to include N A = 0, 1, 2, . . . free copies of protein A and N r = 0, 1 
copies of the same protein, which are bound to the promoter. The master equation ac- 
counts for the temporal variations in the probability distribution. For the autorepressor, 
it takes the form [P2] 

P(N A ,N r ) = gS Nrfi [P(N A -l,N r )-P(N A ,N r )} 

+ d[{N A + 1)P(N A + 1, N r ) - N A P{N A , N r )} 

+ a [S Nr A(N A + l)P(N A + l,N r - 1) - 5 Nrt0 N A P{N A , N r )] 

(5) + ai[S Nr , P(N A -l,N r + l)-S Nrtl P(N A ,N r )}, 

where the g term accounts for the production of proteins and the d term accounts for 
their degradation. The ao (ai) terms describe the binding (unbinding) of proteins to 
(from) the promoter site. 

In numerical integration, the master equation must be truncated in order to keep 
the number of equations finite. This is done by setting a suitable upper cutoff, 7V™ ax , 
on the population size of the free proteins. In order to maintain the accuracy of the 
calculations, the cutoff should be chosen such that the probability of population sizes 
beyond it will be sufficiently small. The master equation exhibits a single steady state 
solution, characterized by P(N A ,N r ) = for all N A and N r . This solution is always 
stable [57] . The average concentration of free A proteins can be obtained from 

Nr* l 

(6) (N A )= E N A P(N A ,N r ). 

N A =0 N t =0 

Other properties of the distribution, such as the variance, can be obtained from the 
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Figure 2: Results for the autorepressor circuit, (a) The concentration of A proteins 
vs. time, obtained from the rate equations and from Monte Carlo simulations. The rate 
equation results quickly reach a steady state. The Monte Carlo results fluctuate around 
this steady state, (b) The steady state probability distribution P(Na) for a cell to contain 
Na copies of protein A, obtained from the master equation. The parameters used are 
g = 0.05, d = 0.001, a = 0.01, a x = 0.01 and d r = (sec- 1 ). 

calculation of higher moments. 

Another useful approach to the study of stochastic dynamics is provided by Monte 
Carlo methods [3 El [16]. The simulation dynamics is Markovian. At each instant, the 
next process to take place is chosen randomly from all the possible processes, where 
each process is assigned with a suitable weight, proportional to its rate. The elapsed 
time is then updated accordingly. Unlike the master equation, which accounts for the 
entire distribution, Monte Carlo simulations follow the temporal variations in protein 
concentrations in a single cell. 

In Fig. [2][a) we present the concentration of A proteins vs. time, obtained from the 
rate equations [Eq. Q] and from Monte Carlo simulations. In Fig. [2jb) we show the 
probability distribution P(Na obtained from the master equation under steady state 
conditions. 

3. The Genetic Switch The genetic toggle switch consists of two proteins, A and B, 
which negatively regulate each other at the transcriptional level [Fig.QJb)]. This architec- 
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ture may lead to two steady states, one dominated by A proteins and the other dominated 
by B proteins. When the population of A proteins is much larger than that of B proteins, 
the A proteins tend to suppress the production of B proteins. Under these conditions, the 
production of A proteins is enhanced, because the declining concentration of B proteins 
is not sufficient to suppress it. Therefore, the system approaches a state rich in A proteins 
and poor in B proteins. Similarly, the system may approach a state rich in B proteins 
and poor in A proteins. Transitions between the two states may take place in response to 
a suitable external signal. Spontaneous transitions, due to random fluctuations, are also 
possible. To qualify as a switch, the system should exhibit bistability. In the deterministic 
description, bistability is defined as the existence of two stable steady state solutions of 
the rate equations. This description does not account for the possibility of spontaneous 
transitions between the two states. In the stochastic description, spontaneous transitions 
are taken into account. Therefore, the condition for bistability is that the rate of spon- 
taneous transitions (due to random fluctuations rather than an external signal) is much 
lower than the rates of all other relevant processes in the system. 

Genetic switch systems exist in nature, and give rise to different cell behaviors in 
different situations. A notable example is the phage A switch 24 . This switch appears 
in A phages, which infect E. coli bacteria and can exist in two exclusive states, one 
called lysogeny and the other called lysis. In addition to natural switches, a synthetic 
switch was constructed and studied in E. coli [BJ. Numerous studies, using rate equations, 
have concluded that cooperative binding is a necessary condition for the emergence of 
bistability 3, 6, 28, 29, 30 . Stochastic analysis reveals the reason to this fact. For a switch 
without cooperative binding, three peaks are obtained in the probability distribution 
function. These peaks corresponds to three possible states for the system: one in which 
A is highly expressed, a second in which B is highly expressed and a third in which both 
proteins are suppressed (a 'deadlock' situation) [13l [14] . Monte Carlo simulations show 
rapid transitions between these three states. The possibility of simultaneous suppression 
of both proteins, prevents the system from functioning as a switch. It causes the system 
to exhibit three states of limited stability instead of the two stable states that are desired. 

It is found that in switch systems in which the A and B repressors exhibit cooperative 
binding, the deadlock situation is removed and bistability emerges. This can be explained 
as follows. The deadlock situation results from a simultaneous binding of A and B repres- 
sors to the corresponding promoter sites. Without cooperative binding, it is sufficient for 
the minority protein to recruit a single copy that will bind and suppress the production 
of the dominant protein. In the case of cooperative binding (for example, with n = 2) 
the minority protein needs to recruit two copies that will bind simultaneously in order 
to suppress the production of the dominant protein. This is much less likely. Therefore, 
cooperative binding induces bistability and enables the system to function as a switch. 

Apart from cooperative binding, there are several other mechanisms that may stabilize 
bistability in genetic switch systems. The most obvious is the exclusive switch, where 
there is an overlap between the promoters of A and B, leaving no room for both to be 
occupied simultaneously. Such situations are encountered in nature, for example, in the 
lysis- lysogeny switch of phage A [23]. It was shown that in the presence of cooperative 
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Figure 3: The probability distribution P(Na,N&) of the concentrations of A and B 
proteins, for the switch with degradation of bound repressors, obtained from the master 
equation. Two sharp peaks are observed, one dominated by A proteins and the other 
dominated by B proteins. The peaks are sharp and are separated by a region with van- 
ishing probabilities. As a result, the transition rate between the two peaks is low and the 
switch is stable. 

binding, the exclusive binding of A and B repressors enhances the stability of the genetic 
switch [29j . This is because in exclusive binding the access of the minority specie to 
the promoter site is blocked by the dominant specie. For the exclusive switch without 
cooperative binding, rate equations exhibit only a single steady state solution, namely, 
the system is not a switch in the deterministic framework. However, in the stochastic 
framework the system exhibits bistability and functions as a switch. The distribution 
P(Na, Nb) exhibits two peaks, one dominated by A proteins and the other dominated 
by B proteins. The exclusive binding prevents the possibility of a deadlock situation. 

In addition to the exclusive switch, there exist other variants of the genetic switch 
circuit (we focus here on systems without cooperative binding). Consider a switch in 
which not only free proteins, but also bound proteins experience degradation. Bound- 
repressor degradation tends to prevent the deadlock situation in which both A and B 
repressors are bound simultaneously. This is due to the fact that degradation removes the 
bound repressor from the system, unlike unbinding, where the resulting free repressor may 
quickly bind again. It turns out that degradation of bound repressors induces bistability 
not only in the stochastic framework but also at the level of deterministic rate equations. 
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Figure 4: The concentration, Na, of A proteins vs. time, obtained from Monte Carlo 
simulations for the switch with degradation of bound repressors. The two states are clearly 
observed: one in which Na fluctuates around 50, and another in which it is nearly zero. 
Transitions between these states occur at an average rate of one transition every ~ 10 6 
sec (about 10 days). 



In Fig. [3] we present the probability distribution P(N\, N&) of the concentrations 
of A and B proteins, for the switch with degradation of bound repressors. Two sharp 
peaks, well separated from each other are observed, illuminating the bistable nature of 
the systems. The parameters used in Fig. [3] and in the rest of the paper are: g = 0.15, 
d = 0.003, a = 0.5, a x = 0.01 and d r = 0.003 (sec- 1 ). 

In Fig. [4] we present the temporal variations of the concentration, Na, of A proteins, 
obtained from Monte Carlo simulations for the switch with degradation of bound repres- 
sors. Two states are clearly observed: one in which Na is dominant and another state 
in which it is suppressed. Note that in spite of the very large fluctuations, the switch is 
stable and the average time between spontaneous transitions is about 10 days. 

Another variant of the genetic switch, which exhibits bistability even at the level of 
rate equations, involves protein-protein interactions, where A and B proteins bind to 
each other and form a complex that does not function as a transcription factor. This 
additional process contributes to the stability of the switch because in such 'pair anni- 
hilation' processes the minority protein is affected more strongly. It is thus less likely to 
bind and suppress the production of the dominant protein. 

4. The Repressilator The repressilator circuit consists of three transcription factors, 
A, B and C, which negatively regulate each other's synthesis in a cyclic manner [Fig.QJc)] . 
This circuit was synthetically constructed on plasmids in E. coli and was found to exhibit 
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Figure 5: The concentration, Na, of A proteins vs. time, obtained from Monte Carlo sim- 
ulations of the repressilator circuit. The oscillations are noisy. Their period and amplitude 
vary from cycle to cycle. 

oscillations in the concentrations of the three transcription factors. To understand the 
origin of these oscillations, consider a situation in which the number of A proteins is 
large. In this case it is likely that one of the A proteins will bind the to b promoter and 
will repress the production of B proteins. The reduced level of B proteins will enable the 
gene c to be fully expressed and the number of C proteins will increase and will start to 
repress gene a. As a result, the number of A proteins will decrease, and gene b will be 
activated, completing a full cycle. The order of appearance of the dominant protein type 
in this cycle is A — > C — > B — > A. 

From a theoretical point of view, oscillations in this system can be obtained (under 
some conditions) both in rate equations and in Monte Carlo simulations [Fig. |5]. The 
oscillations obtained from the rate equations arc regular, forming a stable limit-cycle. The 
oscillations obtained from the Monte Carlo simulations are noisy and irregular. Moreover, 
the period and amplitude differ significantly between the rate equations and the Monte 
Carlo simulations [15] . 

The repressilator system was constructed synthetically on plasmids. When the number 
of plasmids in a cell is small, the numbers of promoter sites and bound transcription 
factors are also small. As a result, one expects large fluctuations in the transcription 
rates regulated by these transcription factors. In this case stochastic methods are required. 
However, when the number of plasmids is large, fluctuations are reduced and the rate 
equations become applicable. Therefore, by gradually changing the number of plasmids in 
the cell one can explore the cross-over from the stochastic regime, where the oscillations 
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are noisy, to the deterministic regime, where the oscillations are regular. 

5. Summary We have presented deterministic and stochastic analysis of three simple 
genetic circuits which involve transcriptional regulation and feedback: the autorepressor, 
the switch and the repressilator. Such systems can be simulated using rate equations, that 
account for the concentrations of the mRNAs and proteins produced by these genes. Rate 
equations are suitable when these concentrations are large and fluctuations are negligible. 
However, when some of the transcription factors and their binding sites in a cell appear in 
low copy numbers, fluctuations become important and the rate equations fail. In this case 
stochastic methods such as the master equation or Monte Carlo simulations are required. 
We have shown that fluctuations give rise to quantitative and qualitative changes in 
the dynamics of the systems. In particular, we demonstrated the fluctuations-induced 
bistability in the exclusive switch and the noisy oscillations obtained in the repressilator 
circuit. 
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